#ifndef _CIRCLE_SYSTEM_H_
#define	_CIRCLE_SYSTEM_H_

#include <algorithm>
#include <cmath>
#include <vector>

#include "common/particle_circle.h"
#include "common/potential_energy.h"
#include "common/circle_union.h"
#include "common/circle_diagram.h"


namespace Packing_common
{

	class Circle_system : public System_potential_energy
	{
	protected:
		typedef Particle_circle::Point_2 Point_2;
		typedef Particle_circle::Vector_2 Vector_2;

	public:

		Circle_system(const std::vector<Particle_circle>& particles, const Container_box& bndbox);

		virtual int optimize(double *trans, bool apply_constraints = false);
		//int circle_optimize(double* trans);
		// objective functions
		//double system_energy_at(const double *const x);
		virtual double gravity_energy_at(const double *const x);
		virtual double molecular_energy_at(const double *const x);

		virtual ~Circle_system() { }

	private:
		// constraint functions
		virtual void constraint_at(const double *const x, double *constraint_values);
		virtual void constraint_grad_at(const double *const x, double *constraint_grad);
		virtual void energy_gradient_at(const double *const x, double *gradient);

	protected:
		const std::vector<Particle_circle>& particle_sys; // initial system state
	};

	class Circle_union_system : public System_potential_energy
	{
		typedef Circle_union::Point_2	Point_2;
		typedef Circle_union::Vector_2	Vector_2;

	public:
		Circle_union_system(const std::vector<Circle_union>& particle_unions, const Container_box& bndbox);

		//Circle_union_system(const std::vector<Circle_union>& particle_unions, const Container_box& bndbox, std::vector<Packing_polygon*>* polygons);

		int optimize(double *trans, bool apply_constraints = false);

		//int circle_union_optimize(double *trans);
		virtual double gravity_energy_at(const double *const x);
		virtual double molecular_energy_at(const double *const x);
		double vector_align_energy(const double *const x);

		double system_energy_at(const double *const x);
		//{
		//	return gravity_energy_at(x) + molecular_energy_at(x) + vector_align_energy(x);
		//}
		~Circle_union_system() { }
		// objective function
		//double system_energy_at_translation(const double *const x);
		//double gravity_energy_at_translation(const double *const x);
	private:
		// constraint functions
		void constraint_at(const double *const x, double *constraint_values);
		void constraint_grad_at(const double *const x, double *constraint_grad);
		void energy_gradient_at(const double *const x, double *gradient);

		void update_circles(const double *const x);

		static int update_interim(const int evalRequestCode, const int n, const int m, const int nnzJ, 
							const int nnzH,	const double * const x,	const double * const lambda, 
							double * const obj, double * const c, double * const objGrad, double * const jac,
							double * const hessian,	double * const hessVector, void * userParams) ;
	private:
		const std::vector<Circle_union>& particle_union_sys;
		std::vector< std::pair<unsigned int, unsigned int> > nbpr;
		double dc;
		double vec_penalty_weight;
		std::vector< std::vector<Particle_circle> > updated_circles;

		//debug
		//std::vector<Packing_polygon*>* pp;
		std::ofstream ofs;
	};
}

#endif